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INTRODUCTION 


In  recent  years  finite  element  method  has  been  widely  used  for 

solving  complex  nonlinear  problems  and  many  large-scale  general  purpose 

computer  programs  have  been  developed.1  The  MARC  and  ANSYS  systems 

have  found  wide  applications,  yet  they  are  quite  expensive.  The  piece- 

wise  linear  analysis  option  of  the  NASTRAN  program  provides  an  algorithm 

2 

for  solving  nonlinear  problems  in  material  plasticity.  The  load  is 
applied  in  increments  such  that  the  stiffness  properties  can  be  assumed 
to  be  constant  over  each  increment.  However,  the  usefulness  of  this 
option  is  quite  limited  because  only  a  few  elements  have  been  imple¬ 
mented.  These  include  rod,  tube,  bar  elements  for  one-dimensional 
problems  and  plate  elements  for  two-dimensional  plane  stress  problems. 
This  paper  describes  the  implementation  of  a  trapezoidal  ring  element 
in  NASTRAN  for  solving  elastic-plastic  problems  of  rotational  symmetry. 

The  theoretical  basis  of  our  implementation  follows  the  approach 

3 

first  developed  by  Swedlow.  A  unique  relationship  between  the 
octahedral  strain  is  assumed  to  exist.  The  material  is  assumed  to 
obey  the  Mises*  yield  criterion,  isotropic  hardening  rule  and  the 
Prandtl-Reuss  flow  relations.  The  explicit  expressions  for  axisym- 
metric  plasticity  are  derived.  The  element  stiffness  matrix  and  the 

■*-W.  Pilkey,  K.  Saczalski,  and  H.  Schaeffer,  "Structural  Mechanics 
Computer  Programs,  Surveys,  Assessments  and  Availability,"  University 
Press  of  Virginia,  1974. 

2"NASTRAN  Theoretical  Manual,"  NASA  SP-221  (01),  1972. 

^  J .  L.  Swedlow,  "The  Thickness  Effect  and  Plastic  Flow  in  Cracked 
Plates,"  ARL  65-216,  Aerospace  Research  Laboratory,  Office  of  Aero¬ 
space  Research,"  October  1965. 
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stress  data  recovery  routines  for  the  trapezoidal  ring  element  are 
developed.  Seven  new  subroutines  are  implemented  and  added  to  the 
NASTRAN  code. 

For  the  purpose  of  demonstration,  two  elastic-plastic  test  prob¬ 
lems  are  solved.  The  first  is  an  infinitely  long  tube  under  uniform 

internal  pressure.  The  NASTRAN  results  are  in  excellent  agreement  with 

4 

an  exact  solution  based  on  a  finite-difference  approach.  The  second 
problem  is  a  thick-walled  cylinder  of  finite  length  loaded  over  part 
of  its  inner  surface.  The  NASTRAN  results  are  compared  with  those 
obtained  by  a  two-dimensional  code  with  the  use  of  quadrilateral  ring 
elements.**  A  good  agreement  between  the  two  results  has  also  been 
achieved. 


CONSTITUTIVE  RELATIONS 

3 

Following  the  development  by  Swedlow,  the  constitutive  relations 
to  be  used  in  our  formulation  for  solving  elastic-plastic  problems  of 
rotational  symmetry  will  be  presented  here.  In  the  development,  a 
unique  relationship  between  the  octahedral  stress  and  the  plastic 
octahedral  strain  is  assumed  to  exist  and  the  use  of  ideally  plastic 
materials  is  excluded.  The  total  strain  components  (er,  £q,  ez,  and 
Yrz)  are  composed  of  the  elastic,  recoverable  deformations  and  the 


•^J.  L.  Swedlow,  "The  Thickness  Effect  and  Plastic  Flow  in  Cracked 
Plates,"  ARL  65-216,  Aerospace  Research  Laboratory,  Office  of  Aero¬ 
space  Research,  October  1965. 

4p.  C.  T.  Chen,  "A  Finite-Difference  Approach  to  Axisymmetric  Plane 
Strain  Problems  Beyond  the  Elastic  Limit,"  To  appear  in  Trans.  25th 
Conference  of  Army  Mathematicians,  1979. 

5P.  C.  T.  Chen,  "Numerical  Solution  of  Gun  Tube  Problems  in  the  Elastic- 
Plastic  Range,"  Proceedings  of  the  1977  Army  Numerical  Analysis  and 
Computer  Conference,  pp.  423-439. 
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plastic  portions  (eP,  ep,  eP,  and  y13  ).  The  rates  of  plastic  flow, 

r  9  z  rz 

(eP,  etc.),  are  independent  of  a  time  scale  and  are  simply  used  for 

convenience  instead  of  the  incremental  values.  The  definitions  of  the 

octahedral  stress  T  and  the  octahedral  plastic  strain  rate  eP  in  the 

°  o 

case  of  rotational  symmetry  are 

T0  -  (l/])[(0r-0e)2  *  (°e-az)!  *  (a2-0r)2  *  6xrz2]1/2  (1) 

t?o  -  (1/3)  2  *  (jjj-^)2  *  C^-jP)2  *  (3/2)(yPj2]1/2  (2) 

where  (a  ,  afl,  a  ,  t  )  are  the  nonvanishing  stress  components. 

1  V  Z  1  Z 

A  unique  relationship  between  tq  and  is  assumed  and  there 
exists  a  function,  M^Ct^),  such  that 

2NW  =  V®S  •  (3) 

The  plastic  modulus,  can  be  related  to  the  slope,  E  ,  of  the 

effective  stress-strain  (cr-e)  curve  by 

1  _  _1_  1  fA. 

3Hr(To)  "  Et  '  E  l4J 

where  E  is  the  elastic  modulus  and 

•  o 

Et  =  a/l 

a  =  (3//2)tq 

e  =  Jl  e0  +  ff/E  .  (5) 

The  material  is  assumed  to  obey  the  Mises*  yield  criterion  and  the 
Prandtl-Reuss  flow  rule.  The  matrix  relationship  for  the  plastic  flow 
in  the  case  of  rotational  symmetry  is 
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where  v  is  the  Poisson's  ratio. 


(6) 


(7) 


Sr  “  (2Waz)/3  ’ 

se  =  (2War)/3  » 

Sz  =  (2az-ar-ae)/3  .  (8) 

For  strain-hardening  materials,  Hp  (or  E-jO  t  0»  we  can  obtain  the 
inverse  of  [C]  numerically  and  this  procedure  is  chosen  in  developing 
NASTRAN  program.  For  ideally-plastic  materials,  MT  =  ET  =  0,  matrix 
[C]  does  not  exist  and  the  NASTRAN  program  fails.  However,  its 
inverse  [C]-*  still  exists  and  the  closed  form  has  been  derived  in 
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Reference  6.  In  the  axisymmetric  case,  the  explicit  form  as  seen  in 
Reference  7  is: 
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where 

A  =  3t02(1  +  Mp/G) 


(9) 


and 

G  =  j  E/(l+v)  .  (10) 

If  we  want  to  remove  the  limitation  that  the  use  of  ideal ly-plastic 
materials  is  excluded,  we  have  to  use  Eq.  (9)  instead  of  Eq.  (7). 


TRAPEZOIDAL  RING  ELEMENT 

The  incremental  displacement  field  employed  for  the  trapezoidal 
ring  element  are 

Au(r,z)  =  3j  +  B2r  +  83Z  +  B4rz  , 

Aw(r,z)  =  B5  +  B6r  +  B?z  +  Bgrz  .  (11) 

The  transformation  from  grid  point  coordinates  to  generalized  coor¬ 
dinates  is 

(B)  =  [r^q]{Aq}  (12) 


.  Yamada,  N.  Yoshimura,  and  T.  Sakurai,  "Plastic  Stress-Strain  Matrix 
and  Its  Application  for  the  Solution  of  Elastic-Plastic  Problems  by 
the  Finite  Element  Method,"  Int.  J.  Mech.  Sci„,  Vol.  10,  pp.  345-354, 
1968. 

n 

P*  C.  T.  Chen,  "Elastic-Plastic  Solution  of  a  Two-Dimensional  Tube 
Problem  by  the  Finite  Element,"  Trans.  Nineteenth  Conf.  of  Army  Mathe¬ 
maticians,  pp.  763-784. 
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where 


{Aq}  =  [Au^  Aw^  Au2,  Aw2>  Au3,  Aw3,  Au4>  Av?4]  , 

(6}T  =  [Blf  e2,  B3,  34,  B5,  36,  B-,  Bg]  .  (13) 

and  the  elements  of  the  inverse  of  the  transformation  matrix  [r„  ]  * 

Bq 

are  the  coefficients  of  the  B's  in  Equations  (11),  evaluated  at  the 
comers  of  the  element. 

The  stiffness  matrix  is  formed  in  the  same  manner  as  that  for  the 
anisotropic  elastic  element.  The  final  form  referred  to  grid  point 
coordinates  is 

w  ■  [WT™'V  ■ 

where 

[K]  =  2-ir/r [B]T [D]  [B]dzdr  .  (14) 

[D] ,  the  matrix  of  material  coefficients,  is  defined  by  Equation  (9) . 
The  [B]  matrix  is  the  same  as  the  elastic  case,  but  now  it  expresses 
the  incremental  strains  in  terms  of  generalized  coordinates 

{Ae}  =  [B]{3>  .  (15) 


NASTRAN  IMPLEMENTATION 

The  NASTRAN  implementation  for  an  elastic-plastic  trapezoidal 
ring  element  follows  the  steps  outlined  in  Section  6.8,  "Adding  a 
structural  element,"  of  Reference  8.  Changes  were  required  in  the 
functional  modules  PLA1,  PLA3,  and  PLA4,  which  included  the  writing 

^"NASTRAN  Programmers'  Manual,"  NASA  SP-223  (01),  1972. 
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of  several  new  subroutines.  These  new  routines  could  easily  be  modeled 
after  the  existing  code  for  the  linear  portion  of  the  program.  There 
are  two  major  differences  in  the  nonlinear  subroutines.  First  is  the 
new  code  for  the  calculation  of  the  material  stiffness  matrix  and 
second  is  that  thermal  stresses  and  element  force  calculations  are 
eliminated  in  the  nonlinear  code. 

The  changes  in  PLA1  allows  this  module  to  identify  the  new  element 
as  a  member  of  the  piecewise  linear  element  set  and  properly  initialize 
the  nonlinear  Element  Summary  and  Element  Connection  Property  Tables. 
Three  element  stress  recovery  subroutines  were  added  to  PLA3:  PSAPRG, 
a  driver  for  stress  data  recovery;  PLTRR1  and  PSTRR2,  phase  I  and  II 
stress  recovery  routines.  Element  stiffness  calculations  in  PLA4 
require  four  new  subroutines:  PKAPRG,  a  driver  for  nonlinear  trape¬ 
zoidal  ring  elements  in  PLA4,  PKTRR1  and  PKTRR2,  stress  recovery 
routines  which  generate  stresses  for  the  computation  of  the  nonlinear 
material  matrix;  and  PKTRAP,  the  stiffness  matrix  generation  routine 
for  nonlinear  trapezoidal  ring  elements. 

The  computer  system  available  for  this  work  is  IBM  360  Model  44. 

In  order  to  add  the  new  code  into  Link  13,  it  became  necessary  to  add 
a  new  branch  to  the  overlay  trees  to  contain  the  new  elements. 
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NUMERICAL  EXAMPLES 


For  the  purpose  of  demonstration,  two  elastic-plastic  problems  of 

rotational  symmetry  were  solved  and  compared  with  other  results.4’5 

The  first  is  an  infinitely  long  tube  under  uniform  internal  pressure 

and  the  second,  a  thick-walled  cylinder  of  finite  length  loaded  over 

part  of  its  inner  surface.  The  material  properties  for  both  problems 

are  the  same.  The  elastic  constants  are  E  =  30xl06  psi,  v  =  0.3  and 

the  effective  stress-strain  curve  is  represented  by  three  line  segments 

connecting  the  four  points  in  the  (e-a)  plane,  (E,a)  =  (0.0, 0.0), 

(0.005,  150,000  psi) ,  (0.055,  225,000  psi),  (0.1,  225,000  psi). 

EXAMPLE  1.  Consider  an  infinitely  long  tube  subjected  to  uniform 

internal  pressure  p.  The  plane  strain  condition  is  assumed.  The  tube 

of  outside  radius  2"  and  inside  radius  1"  has  been  divided  into  25 

trapezoidal  ring  elements.  The  numerical  results  based  on  the  NASTRAN 

program  have  been  obtained.  For  this  problem,  a  new  finite-difference 
4 

approach  can  be  used  to  generate  exact  solution  and  to  assess  the 
accuracy  of  the  NASTRAN  code.  Some  of  the  results  for  the  displace¬ 
ments  and  stresses  are  presented  graphically  in  Figures  1  and  2. 
Twenty-five  load  increments  are  used  in  NASTRAN  as  shown  in  Figure  1. 

The  radial  displacements  at  the  inside  as  well  as  outside  surface  are 

4P.  C.  T.,  "A  Finite-Difference  Approach  to  Axisymmetric  Plane  Strain 

Problems  Beyond  the  Elastic  Limit,"  To  appear  in  Trans.  25th  Confer¬ 
ence  of  Army  Mathematicians,  1979. 

^P.  C.  T.  Chen,  "Numerical  Solution  of  Gun  Tube  Problems  in  the  Elastic- 
Plastic  Range,"  Proceedings  of  the  1977  Army  Numerical  Analysis  and 
Computer  Conference,  pp.  423-439. 
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shown  as  functions  of  internal  pressure.  Figure  2  shows  the  distri¬ 
butions  of  radial,  tangential  and  axial  stress  components  in  a  pressur¬ 
ized  tube  when  half  of  the  tube  is  plastic.  The  pressure  required  to 

achieve  this  state  is  0.7378  an  based  on  NASTRAN  code  and  0.7356  a 

o  o 

4 

based  on  the  finite-difference  solution.  Both  codes  indicate  that 

the  maximum  tangential  stress  occurs  at  the  elastic-plastic  interface. 

As  demonstrated  in  Figures  1  and  2,  the  NASTRAN  results  are  in  excel- 

4 

lent  agreement  with  those  based  on  the  finite-difference  approach. 

EXAMPLE  2.  Consider  a  two-dimensional  elastic-plastic  thick- 
walled  cylinder  problem  as  shown  in  Figure  3.  The  tube  with  inner 
radius  (lff),  outer  radius  (2M)  and  length  (4”)  is  loaded  uniformly 
over  a  middle  portion  (2")  of  the  inner  surface.  The  mesh  generation 
and  the  loading  for  the  half  of  the  undeformed  structure  is  shown  in 
the  figure,,  This  problem  was  solved  in  Reference  5  based  on  a  scale 
loading  approach.  The  first  load  factor  is  the  upper  limit  of  the 
elastic  solution  and  ten  additional  increments  were  needed  until  one 
of  the  outside  elements  becomes  yielded.  The  same  load  factors  were 
used  to  obtain  the  NASTRAN  solution.  Both  programs  indicate  that  the 
sequence  in  which  the  elements  become  plastic  is  1,5,9,2,13,6,10,3, 
7,14,11,17,4.  Some  steps  will  cause  more  than  one  element  to  become 

^P.  C.  T.  Chen,  "A  Finite-Difference  Approach  to  Axisymmetric  Plane 
Strain  Problems  Beyond  the  Elastic  Limit, 11  To  appear  in  Trans.  25th 
Conference  of  Army  Mathematicians,  19790 

5p.  C.  T.  Chen,  Numerical  Solution  of  Gun  Tube  Problems  in  the  Elastic- 
Plastic  Range,"  Proceedings  of  the  1977  Army  Numerical  Analysis  and 
Computer  Conference,  pp.  423-439. 
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plastic  and  those  elements  with  effective  stress  o  _>  0.99  oq  have  been 
considered  as  yielded.  The  numerical  results  for  the  radial  displace¬ 
ment  at  the  inside,  ua  (point  1)  and  outside,  u^  (point  5)  as  functions 
of  internal  pressure  are  shown  in  Figure  4.  The  stress  components  at 
the  centroid  of  one  inside  element  (No.  1)  are  shown  in  Figure  5.  The 
effect  of  loading  history  on  the  displacements  and  stresses  can  be  seen 
from  Figures  4  and  5.  A  comparison  of  the  results  between  NASTRAN 
program  and  reference  5  indicates  that  very  good  agreement  has  been 
achieved. 

CONCLUSION 

An  elastic-plastic  trapezoidal  ring  element  has  been  implemented 
in  NASTRAN  computer  program.  Its  application  to  elastic-plastic  prob¬ 
lems  of  rotational  symmetry  has  been  demonstrated  by  solving  two  thick- 
walled  tube  problems.  The  NASTRAN  results  for  both  problems  are  in 
excellent  agreement  with  the  other  results. 


^P.  C.  T.  Chen,  "Numerical  Solution  of  Gun  Tube  Problems  in  the  Elastic- 
Plastic  Range,"  Proceedings  of  the  1977  Army  Numerical  Analysis  and 
Computer  Conference,  pp.  423-439. 
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Figure  1.  Radial  displacements  in  a  pressurized  tube. 
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Figure  3.  Two-dimensional  elastic-plastic  tube. 
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Figure  4.  Radial  displacements  in  a  two-dimensional  tube. 
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